7-6 shortness of breath

Shortness of breath.

Authors
Affiliations

James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Published

July 26, 2023

Load packages
library(ASCOTr)
library(tidyverse)
library(lubridate)
library(kableExtra)
library(patchwork)
library(cmdstanr)
library(posterior)
library(bayesplot)
library(ggdist)
library(lme4)
library(broom)
library(broom.mixed)
library(bayestestR)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))

bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)
Prepare datasets
all_dat <- readRDS(file.path(ASCOT_DATA, "all_data.rds"))

# FAS-ITT
fas_itt_dat <- ASCOTr:::make_fas_itt_set(all_dat)
fas_itt_nona_dat <- fas_itt_dat |>
  filter(!is.na(out_sob))

# ACS-ITT
acs_itt_dat <- ASCOTr:::make_acs_itt_set(all_dat)
acs_itt_nona_dat <- acs_itt_dat |>
  filter(!is.na(out_sob))

# AVS-ITT
avs_itt_dat <- ASCOTr:::make_avs_itt_set(all_dat)
avs_itt_nona_dat <- avs_itt_dat |>
  filter(!is.na(out_sob))
Load models
logistic_mod <- compile_cmdstanr_mod(
  file.path("binary", "logistic"), dir = "stan")
logistic_site_epoch <- compile_cmdstanr_mod(
  file.path("binary", "logistic_site_epoch"), dir = "stan")
Code
make_summary_table_anticoagulation <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3", "C4"),
      labels = str_replace(intervention_labels()$CAssignment[-1], "<br>", " "))) %>%
    group_by(CAssignment) %>%
    summarise(
      Patients = n(),
      Known = sprintf(
        "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
      Missing = sprintf(
        "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
      `Shortness of breath day 28` = sprintf(
        "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Patients = n(),
          Known = sprintf(
            "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
          Missing = sprintf(
            "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
          `Shortness of breath day 28` = sprintf(
            "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(key, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(font_size = 10, latex_options = "HOLD_position")  %>%
    row_spec(nrow(tab), bold = TRUE)
}


make_summary_table_antiviral <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(AAssignment = factor(
      AAssignment, 
      levels = c("A1", "A2"),
      labels = str_replace(intervention_labels()$AAssignment[-1], "<br>", " "))) %>%
    group_by(AAssignment) %>%
    summarise(
      Patients = n(),
      Known = sprintf(
        "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
      Missing = sprintf(
        "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
      `Shortness of breath day 28` = sprintf(
        "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(AAssignment = "Overall") %>%
        summarise(
          Patients = n(),
          Known = sprintf(
            "%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
          Missing = sprintf(
            "%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
          `Shortness of breath day 28` = sprintf(
            "%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
      )
    ) %>%
    mutate(AAssignment = fct_inorder(AAssignment)) %>%
    gather(key, value, -AAssignment, factor_key = T) %>%
    spread(key, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(font_size = 10, latex_options = "HOLD_position")  %>%
    row_spec(nrow(tab), bold = TRUE)
}



make_out_sob_table_C <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = paste0("C", 0:4),
      labels = intervention_labels()$CAssignment)
    ) %>%
    group_by(CAssignment) %>%
    summarise(
      Randomised = n(),
      `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
      `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
      `Died within 28 days` = sprintf("%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Randomised = n(),
          `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(out_sob)), 100 * mean(is.na(out_sob))),
          `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(out_sob)), 100 * mean(!is.na(out_sob))),
          `Died within 28 days` = sprintf("%i (%.1f)", sum(out_sob, na.rm = TRUE), 100 * mean(out_sob, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(CAssignment, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(
      font_size = 9, 
      latex_options = "HOLD_position")  
}

Descriptive

Relationship anti-coagulation to outcome
save_tex_table(make_summary_table_anticoagulation(
      acs_itt_dat,
      "latex"), 
      "outcomes/secondary/7-6-anticoagulation-summary")
make_summary_table_anticoagulation(
    acs_itt_dat
)
n (\%) Patients Known Missing Shortness of breath day 28
Low dose 610 577 (94.6) 33 (5.4) 115 (19.9)
Intermediate dose 613 584 (95.3) 29 (4.7) 110 (18.8)
Low dose with aspirin 283 271 (95.8) 12 (4.2) 59 (21.8)
Therapeutic dose 50 44 (88.0) 6 (12.0) 11 (25.0)
Overall 1556 1476 (94.9) 80 (5.1) 295 (20.0)
Table 1: Primary outcome by anti-coagulation intervention.
Relationship anti-viral to outcome
save_tex_table(make_summary_table_antiviral(
     avs_itt_dat,
      "latex"), 
      "outcomes/secondary/7-6-antiviral-summary")
make_summary_table_antiviral(
     avs_itt_dat
)
n (\%) Patients Known Missing Shortness of breath day 28
Standard of care 73 67 (91.8) 6 (8.2) 35 (52.2)
Nafamostat 82 73 (89.0) 9 (11.0) 36 (49.3)
Overall 155 140 (90.3) 15 (9.7) 71 (50.7)
Table 2: Primary outcome by anti-viral intervention.

Age

Relationship age to shortness of breath
agedat <- fas_itt_dat %>%
  dplyr::count(out_sob, AgeAtEntry) %>% 
  spread(out_sob, n, fill = 0) %>% 
  mutate(
    n = `0` + `1` + `<NA>`,
    p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p1 <- ggplot(agedat, aes(AgeAtEntry, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry")
p2 <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_vline(xintercept = 60, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\nshortness of breath\nday 28", x = "Age at entry")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-age.pdf"), p, height = 2, width = 6)
p

Figure 1: Relationship (logistic regression linear in age) between age at entry and shortness of breath at day 28.

Sex

Relationship sex to outcome
tdat <- fas_itt_dat %>%
  dplyr::count(Sex, out_sob) %>%
  group_by(Sex) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Sex, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Sex")
p2 <- ggplot(tdat, aes(Sex, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Sex")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-sex.pdf"), p, height = 2, width = 6)
p

Figure 2: Shortness of breath by sex.

Oxygen

Relationship oxygen to outcome
tdat <- fas_itt_dat %>%
  dplyr::count(supp_oxy2 = factor(supp_oxy2, labels = c("Not required", "Required")), out_sob) %>%
  group_by(supp_oxy2) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(supp_oxy2, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Supplemental oxygen")
p2 <- ggplot(tdat, aes(supp_oxy2, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Supplemental oxygen")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-oxygen.pdf"), p, height = 2, width = 6)
p

Figure 3: Shortness of breath by oxygen.

Country

Relationship country to outcome
tdat <- fas_itt_dat %>%
  dplyr::count(Country = factor(
    Country, 
    levels = c("IN", "AU", "NP", "NZ"),
    labels = c("India", "Australia", "Nepal", "New\nZealand")), out_sob) %>%
  group_by(Country) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(tdat, aes(Country, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Country of enrolment")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-country.pdf"), p, height = 2, width = 6)
p

Figure 4: Shortness of breath by country.

Site

Relationship site to outcome
tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(
    Country_lab = Country,
    Site_lab = fct_infreq(Location),
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")),
    Site = PT_LocationName,
    out_sob) %>%
  group_by(Country, Site) %>%
  spread(out_sob, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()
p1 <- ggplot(tdat, aes(Site_lab, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(tdat, aes(Site_lab, p_1)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Site of enrolment") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p <- p1 / p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-country-site.pdf"), p, height = 4, width = 6.25)
p

Figure 5: Shortness of breath at day 28 by site within country.

Calendar Time

Relationship calendar date to outcome
caldat <- all_dat %>% 
  filter_fas_itt() %>%
  dplyr::count(out_sob, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(out_sob, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`),
         n = `1` + `0` + `<NA>`)
p1 <- ggplot(caldat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_point() +
    labs(
      y = "Proportion\nshortness of breath\nday 28", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12) +
  ylim(0, NA)
p2 <- ggplot(caldat, aes(mth, n)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p <- p2 | p1
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-6-calendar-time.pdf"), p, height = 2, width = 6)
p

Figure 6: Relationship between calendar date and the primary outcome.

Analyses

FAS-ITT

Code
res <- fit_primary_model(fas_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "out_sob")
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-6-primary-model-fas-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.71 (0.30, 1.56) 0.76 (0.33) 0.81
Intermediate-dose 0.78 (0.55, 1.10) 0.79 (0.14) 0.92
Low-dose with aspirin 1.17 (0.75, 1.82) 1.20 (0.27) 0.24
Therapeutic-dose 0.97 (0.37, 2.39) 1.07 (0.53) 0.53
Ineligible aspirin 1.19 (0.40, 3.19) 1.34 (0.74) 0.37
Age ≥ 60 2.10 (1.50, 2.94) 2.14 (0.37) 0.00
Female 0.99 (0.72, 1.35) 1.00 (0.16) 0.54
Oxygen requirement 1.23 (0.86, 1.75) 1.25 (0.23) 0.14
Australia/New Zealand 2.67 (0.72, 9.63) 3.32 (2.48) 0.07
Nepal 0.48 (0.12, 2.36) 0.66 (0.74) 0.83
Posterior summaries for model parameters (fixed-effects), FAS-ITT.
Odds ratio summary for epoch and site
p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-epoch-site-terms-fas-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.
Code
p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-fas-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7667 0.8095 0.8068 0.8222 0.7507 0.7393 0.8046 0.7920
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["gamma_epoch"])

Posterior Predictive Checks

Code
y_ppc <- res$drws$y_ppc
ppc_dat <- bind_cols(fas_itt_nona_dat, tibble(y_ppc = y_ppc))

grp_ppc <- function(grp) {
  ppc_dat %>%
    group_by(grp = !!grp) %>%
    summarise(
      mean_y = mean(out_sob),
      rvar_mean_y_ppc = rvar_mean(y_ppc)
    )
}

plot_grp_ppc <- function(dat, lab = "") {
  dat %>%
    ggplot(aes(y = grp, xdist = rvar_mean_y_ppc)) +
    stat_interval(size = 2) +
    geom_point(aes(x = mean_y, y = 1:nrow(dat)), colour = "red", shape = 23) +
    labs(
      x = "Posterior predictive\nproportion", 
      y = lab)  
}

ppc_A <- grp_ppc(sym("AAssignment"))
ppc_C <- grp_ppc(sym("CAssignment"))
ppc_ctry <- grp_ppc(sym("Country"))
ppc_epoch <- grp_ppc(sym("epoch"))
ppc_site <- ppc_dat %>%
  group_by(Country = Country, grp = site) %>%
  summarise(
    mean_y = mean(out_sob),
    rvar_mean_y_ppc = rvar_mean(y_ppc)
  )

p0 <- plot_grp_ppc(ppc_A, "Anti-\nviral") + labs(x = "")
p1 <- plot_grp_ppc(ppc_C, "Anti-\ncoagulation") + labs(x = "")
p2 <- plot_grp_ppc(ppc_ctry, "Country") + labs(x = "")
p3 <- plot_grp_ppc(ppc_epoch, "Epoch") + labs(x = "")
p4 <- plot_grp_ppc(ppc_site %>% filter(Country == "IN"), "Sites\nIndia")  + labs(x = "")
p5 <- plot_grp_ppc(ppc_site %>% filter(Country == "AU"), "Sites\nAustralia") + labs(x = "")
p6 <- plot_grp_ppc(ppc_site %>% filter(Country == "NP"), "Sites\nNepal")
p7 <- plot_grp_ppc(ppc_site %>% filter(Country == "NZ"), "Sites\nNZ")
p <- ((p3 | (p0 | p1) / p2) / 
        ( (p4 | p5) / (p6 | p7) + 
            plot_layout(heights = c(3, 1)) ) ) +
  plot_layout(heights = c(1, 1.5), guides = "collect") &
  theme(legend.position = "bottom")
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-6-primary-model-fas-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.5)
p

ACS-ITT

Fit primary model
res <- fit_primary_model(acs_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "out_sob")
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-6-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.85 (0.34, 2.07) 0.94 (0.45) 0.64
Intermediate-dose 0.77 (0.54, 1.09) 0.78 (0.14) 0.93
Low-dose with aspirin 1.17 (0.75, 1.81) 1.20 (0.27) 0.24
Therapeutic-dose 0.95 (0.36, 2.35) 1.06 (0.52) 0.54
Ineligible aspirin 1.23 (0.41, 3.34) 1.39 (0.77) 0.36
Age ≥ 60 2.04 (1.45, 2.90) 2.07 (0.37) 0.00
Female 1.02 (0.75, 1.40) 1.04 (0.16) 0.44
Oxygen requirement 1.28 (0.88, 1.85) 1.30 (0.25) 0.10
Australia/New Zealand 2.28 (0.61, 8.37) 2.84 (2.16) 0.11
Nepal 0.50 (0.12, 2.46) 0.69 (0.69) 0.82
Posterior summaries for model parameters (fixed-effects), FAS-ITT.
Odds ratio summary for epoch and site
p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-epoch-site-terms-acs-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.
Code
p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-acs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7511 0.7915 0.7607 0.7888 0.7786 0.7750 0.7824 0.8094
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["gamma_epoch"])

AVS-ITT

Code
res <- fit_primary_model(
  avs_itt_nona_dat, 
  logistic_mod, 
  ctr = contr.equalprior,
  outcome = "out_sob", 
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile"), 
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5)
)
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Oxygen requirement", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)")
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-6-primary-model-avs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.76 (0.35, 1.62) 0.82 (0.33) 0.76
Intermediate-dose 0.78 (0.32, 1.88) 0.86 (0.41) 0.71
Low-dose with aspirin 0.58 (0.13, 2.38) 0.75 (0.61) 0.77
Therapeutic-dose 1.65 (0.51, 5.57) 2.00 (1.36) 0.20
Age ≥ 60 3.04 (1.33, 7.34) 3.37 (1.57) 0.00
Female 1.16 (0.55, 2.51) 1.26 (0.51) 0.35
Oxygen requirement 2.56 (1.09, 6.20) 2.83 (1.32) 0.01
CRP (2nd tertile) 0.83 (0.32, 2.18) 0.94 (0.49) 0.65
CRP (3rd tertile) 0.73 (0.28, 1.89) 0.82 (0.42) 0.74
CRP (unknown) 1.08 (0.30, 3.88) 1.34 (0.98) 0.45
Posterior summaries for model parameters (fixed-effects), AVS-ITT.
Code
p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-6-primary-model-avs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0006 0.9552 0.9447 0.9590 1.0144 0.9676 0.9890 0.9826
Code
mcmc_trace(res$drws["beta"])

End of script.